tolXI = 1e-11; maxit = 5000;
[V,e] = eig(thch.Pmat');
ie    = find(diag(e) >= 1-1e-12);
XI_ss  = V(:,ie)/sum(V(:,ie)); 
XI_err = 8; it = 1;
while (XI_err > tolXI) && (it <= maxit)    
    XI_new = thch.Pmat'*XI_ss;
    XI_err = norm(XI_new - XI_ss);
    XI_ss  = XI_new;    
    if it == maxit
        disp(['couldnt find X_ss, XI_err = ', num2str(XI_err)])        
    end    
    it = it + 1;    
end